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Abstract 

A coarse-grained particle model for incompressible Navier-Stokes (NS) equation is 
proposed based on spatial filtering by utilizing smoothed particle hydrodynamics 
(SPH) approximations. This model is similar to our previous developed SPH dis- 
cretization of NS equation (Hu X.Y. & N.A. Adams, J. Comput. Physics, 227:264- 
278, 2007 and 228:2082-2091, 2009) and the Lagrangian averaged NS (LANS-a) 
turbulence model. Other than using smoothing approaches, this model obtains par- 
ticle transport velocity by imposing constant a which is associated with the particle 
density, and is called SPH-o" model. Numerical tests on two-dimensional decay and 
forced turbulences with high Reynolds number suggest that the model is able to 
reproduce both the inverse energy cascade and direct enstrophy cascade of the ki- 
netic energy spectrum, the time scaling of enstrophy decay and the non-Guassian 
probability density function (PDF) of particle acceleration. 
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1 Introduction 



The smoothed particle hydrodynamics (SPH) method is a fully Lagrangian, 
grid free method. Since its introduction by Lucy [14j and Gingold and Mon- 
aghan [B], SPH has been applied to a wide range of macroscopic and meso- 
scopic flow problems [I9][l0]. Though the SPH method was originally de- 
veloped for astrophysical problems involving compressible fluids, it has been 
extended to problems involving incompressible fluids by using either a weakly 
compressible model of the fluid [T7], or by algorithms designed to solve the 
full incompressible equations [3] [11] [12]. 

Many of the incompressible flow problems, such as flood and coastal flows, to 
which SPH has been applied are turbulent. Since the direct numerical simula- 
tion of these problems is not always feasible, turbulence modeling is required 
for the computational more efficient coarse-grained numerical simulation. One 
straightforward approach of SPH turbulence modeling is applied the turbu- 
lence models originally developed for Euelrian methods directly [7] [25] [27] . 

Monaghan [18] first noticed the similarity between the version of SPH called 
XSPH [16] and the Lagrangian averaged Navier-Stokes (LANS) turbulence 
model [8] [9] on the relation between the velocity determined from particle 
momentum (momentum velocity) and the transport velocity, and proposed 
a turbulence model specifically for the SPH method. In this model, the SPH 
particle moves with the transport velocity smoothed from momentum velocity 
by an iterative algorithm and a dissipation term is introduced to mimic the 
standard large eddy simulation (LES) model originally developed for Eulerian 
methods. A further modification of this model (SPH-e) is to obtain trans- 
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port velocity directly by the XSPH method with a parameter e pO]. On the 
other hand, we have noticed the importance of SPH particle moving with the 
velocity different from its momentum velocity when simulating flows beyond 
small Reynolds number in our previous developed incompressible SPH method 
pT| [12] . In this method, other than using the XSPH method or smoothing ap- 
proaches, the Eulerian incompressibility condition (zero velocity divergence) 
and the Lagrangian incompressibility condition (constant density) are used 
respectively to determine the momentum velocity and the transport velocity. 

In this paper, we propose a coarse-grained particle system for turbulence sim- 
ulation based on spatial filtering the Navier-Stokes (NS) equation by utilizing 
SPH approximations. Since the resulting particle equations are similar to those 
of the above mentioned incompressible SPH method, except for an additional 
effective stress term introduced by moving particle with transport velocity, 
the same numerical method is applied. The numerical tests show that, while 
achieving good accuracy for resolved flow, the present model can recover the 
spectral and statistical properties of the two-dimensional decay and forced 
turbulences with high Reynolds number. 



2 Model 



We consider the incompressible isothermal NS equation in Lagrangian form 



d\ dv 

11= m^^-^^ 

dp 
'dt 



P 



or V ■ V = 0, 



(1) 
(2) 
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where p, p and v are fluid density, pressure and velocity, respectively, and 
u = rj/p is the kinematic viscosity. Note the two expressions (constant density 
and zero velocity-divergence) in Eq. are formally equivalent. 

The equation of the LANS-a model [H] [S] is a coarse-grained equation, written 
in the form of Eulerian mean or filtered velocity v and Lagrangian mean or 
transport velocity v as 



_ + v.Vv 



1. 



-Vp + i/V^v, 
P 

V ■ v = V ■ V = 0, 



v=(l -a^V^)v. 



(3) 

(4) 
(5) 



The Helmholtz operator (1 — a^V^) in Eq. ([5]) suggests that the transport 
velocity v is smoother than the filtered velocity v. One can think of the pa- 
rameter a as the length scale associated with the width of an extra filter which 
smooths V to obtain v. This has been done explicitly in the SPH-e model [20] , 
where v is based on the filtering 

V = V + e y" v(r') - v(r)G'(r' - r, l)dr', (6) 

where G is the filter with width / and e is a constant parameter. The basic idea 
of smoothing is to prevent the production of small length-scale flow structures. 

In the next section, we propose a coarse-grained particle system based on 
filtering the NS equation with numerical techniques in the SPH method. This 
approach is different from the SPH-e model, which is devised from the SPH 
discretized form of Eckart's Lagrangian [4J. 
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2.1 Coarse-grained NS equation and SPH method 



Assume that the incompressible flow fleld is coarse-grained into a particle 
system with spatial filtering, the variables on particles are obtained by 

^, = Gi*i; = J ijW{r - r„ h)dr, (7) 

where and h are the center and width, respectively, of the filter, W{r — ri,h), 
which is radially symmetric respect to and has the properties J W{r — 
Ti, h)dr = 1, / VW{r - ri)dr = and lim/j^o W{r - r^, h) = 6{r - Fj). Note 
that this filter can also take the role as the SPH kernel function with smoothing 
length h. Substitute with the coordinate r into Eq. ([7]), one has the particle 
position Fj at the center of the filter. The motion of particle is determined by 
its transport velocity, i.e. 

Note that Vj can be different from the filtered or momentum particle velocity 
Vj, which is obtained by substitute ip with velocity v into Eq. ([7j). Similarly, 
the filtered derivatives on particles can be written as 

= - J ipVW{r-Vi,h)dT, (9) 

= ^ - / V-v ■ VW{r - r„ h)dv, (10) 

Since only the filtered values ipi are known, it is impossible to apply the exact 
filtering Gi * Vtp and Gi * dtp/dt, an approximated filtering is carried out after 
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the flow fleld is flrst reconstructed 

where a = X^/c ^^(1" — ffc) is a measure of local length scale, which is larger in 
a dense particle region than in a dilute particle region [13] . From the identity 
J2j W{r — rj)/a = 1, the total volume V can be written as 

^ = E = E / -W{r - r,)dr (12) 

where aj equals to the inverse of particle volume approximately, i.e. cxj = 
Y.kW{Yj - Tfc) = EfcW^jfc ~ With Eq. © and by using Eq. ^ and 

the properties of the kernel function, the approximation of Gi * Vip can be 
obtained by 



G,*Vij^-Y. / -Vl^(r - Vi)W{v - v,)dv ^ -^j^W.j, (13) 

J 3 J 



where ip^^ = {jpi + iljj)/2 is the average between particle i and particle j. 
Note that Eq. ( fT3l) is a typical SPH discretization of the gradient operator. 
Substitute ip with pressure p and velocity v, respectively, in Eq. (fT3ll . the 
filtered pressure gradient and velocity divergence can be approximated as 

G,*Vp^Y. -P^J^WiJ^ * V ■ V ^ 5] -V,, Viy.,. (14) 



Similarly, if the average directional derivative is approximated as VV'jj ~ 
Gijiipi — where and rjj are unit vector and distance, respectively, 

from particle i to particle j, the filtered velocity Laplacian can be approxi- 
mated with Eq. (fT3|) as 

Gi * VV = * V ■ Vv ^ ^ — — ■ VWij. (15) 
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Using Eq. ( !T2|) . the particle density is given by 

Pj = (16) 

where is the constant mass of a particle. The variation of the particle 
density is 

^^^,,J:U,.VW,, (17) 

where Vjj = Vj — is the difference of transport velocity between particle i 
and particle j. Using particle transport velocity Vj, Gi*dil)/dt can be rewritten 



as 



= ^ + Vi ■ + j ^(vi - v) ■ WW{v - r^, h)dv 

« ^ + E v.. ■ VH^.„ (18) 

where ipij = i^i — i^j- Substitute ip with momentum density pv in Eq. (fT8|) . 
one has 

G..^«%i + i:^(H«v«.V,y„, (19) 

where {p'y)ij = PiV^ — pj^j, and the second term is an additional effective stress 
term. With Eqs. ([9]), (fT7|) and (fT9|) . Eqs. ([1]) and ([2]) can be coarse-grained into 



dt ■ ffiffj 



,(20) 



E ■ viy., = E ;^v,, ■ vi^,, = o. (21) 

The particle system of Eqs. (E]), (HZD, (EO]) and (EH) is very similar to the SPH 
discretization of NS equation used in Hu and Adams [11] [12], except that 
there is an additional effective stress term in Eq. (1201) . Actually, without the 
additional effective stress term, Eqs. (1201) and ( |2T1) are a SPH discretization 
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of the Eqs. ([3]) and (jl]) of the LANS-a model. Since this term is induced by 
filtering, its contribution is expected to be small when the flow is well resolved 
as the influence of filtering is negligible. Note that this term is similar to the 
additional effective stress term in the SPH-e model, where only the filtered 
velocities are token into account together with the tunable parameter e. 

To solve this coarse-grained system, the transport velocity Vj is obtained from 
the filtered velocity Vj. For the present particle system, since the local length 
scale is characterized by cTj which variates according to Vj, preventing the 
production of small length scale is equivalent to constrain the constant density 
condition. This is already included in the present system as the first expression 
of Eq. (pTj) . Similar to the denotation of the SPH-e model, in which Vj is 
obtained by a smoothing approach with a parameter e, the present particle 
system is denoted as the SPH-cr model since Vj is obtained by constraining o",. 

3 Numerical method 

As mentioned above that the presented particle system is very similar to the 
SPH discretization of NS equation used in Hu and Adams [TT][T2] . the same 
fractional time-step integration approach can be used. First, the transport 
velocity or half-time-step velocity is obtained by 



where At is time-step size and f is the sum of the viscous term and the ad- 
ditional effective stress term. The pressure field p[ is to impose the constant 
density condition, i.e. the first expression in Eq. (1211) . Subsequently, the par- 




(22) 
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tide position at the new time step is calculated by 

r^' = + ViAt, (23) 
and the particle velocity at the new time step is obtained by 

^r' = < + (fr - ^vp.^ At, (24) 

where the pressure field Pi is to impose the zero velocity-divergence, i.e. the 
second expression in Eq. f l2T]) . 

3.1 Constant- density condition 

We split Eq.f l22p into an intermediate step and a correction step. An interme- 
diate velocity v* is obtained by 

v* = vr + ^f.(r",v"), (25) 

and the intermediate pressure is obtained by solving the discretized Poisson 
equation 

where p* = |p'At^, is the initial value and the left-hand-side is 



(27) 

Note that Eq. fl26l) takes account the accumulated density error, which is given 
by the second term on the right-hand-side. The generalized minimum residual 
(GMRES (g), q is the number of inner iteration steps without restarting) 
method is used to solve Eq. ( 126|) with the convergence condition min(i?ej) < e, 
where Rci is residue, suggesting that the maximum relative variation of (Tj is 



9 



about e, typically from 1% to 0.1%. With Eq. (123|) . the new particle position 
can be obtained by 

= r7 + Atv* - -1 5: -^P- • VW^.,. (28) 

3.2 Zero velocity- divergence condition 

An intermediate velocity at the full time step v*'"+^ is obtained by 

v*'"+^ = VI + Atf, (r", v") . (29) 
The pressure is obtained by solving the discretized Poisson equation 

where p** = pAt. Note that, different from the exact projection for the 
const ant- density constraint, the zero- velocity-divergence constraint is enforced 
by approximate projection, see details in Hu and Adams [11][T2]. The velocity 
at the full time step v"^^ is obtained by 

n+l ^ *,n+l _J_y- J_p-. (31) 

3.3 Time-step criteria 

Besides viscous- diffusion condition, different from Hu and Adams [11] [12] , in 
the present work, the so-called Lipshitz or deformational CFL (DeCFL) con- 
dition is used as the stability condition 

/ ^ 

At < 0.25 min ——i, , — . (32) 

V Vv,: La,, u 
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where ||Vvj||max gives the magnitude of the maximum strain rate of the flow. 
The DeCFL condition suggests that the trajectories of the particles can not 
cross each other. Numerical experiments show that the DeCFL condition al- 
lows larger time-step size than the CFL condition if the flow is moderate. 
However, for flow with high shear rate, the two conditions gives comparable 
time-step sizes. Also note that, though large time-step size according to Eq. 
(1321) does not lead to numerical instability, it increases the temporal truncation 
error. 

4 Numerical simulations 

Several two-dimensional test cases are provided to assess the potential of the 
SPH-o" model. First, the performance of the model for recovering flow with 
moderate Reynolds number is tested. Then two-dimensional decay and driven 
turbulent flows with high Reynolds number are tested. In order to achieve suf- 
ficient high Reynolds number, the physical viscosity in NS equation is switched 
off. For all cases a quintic spline kernel [21] is used as smoothing function. A 
constant smoothing length h, which is kept equal to the initial distance be- 
tween the neighboring particles, is used for all test cases. As elliptic solver a 
diagonally preconditioned GMRES (g) method is used. For all test cases, the 
computation is performed on a domain < a; < 1 and < y < 1 with peri- 
odic boundary conditions in both directions. For spectrum analysis with fast 
Fourier transform (FFT), the values on a uniform grid are reconstruct from 
the particles by a remesh method with a 4th order interpolation kernel [1], 
which is much more computational efficient than the SPH Fourier transform 
used in Robinson and Monaghan [21]. 
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4-1 Two-dimensional Taylor- Green flow 

The two-dimensional viscous Taylor-Green flow is a periodic array of vortices, 
where the velocity 

u = -Ue^^ cos(27rx) sin(27rt/), v = f/e^* sin(27rx) cos(27r?/) (33) 

is an exact solution of the incompressible NS equation, b = — Svr^/Re is the 
decay rate of velocity field. We consider a case with Re = 100, which has 
been used to test different incompressible SPH methods ^ [5] [11] [12] . The 
computational setup is the same as that of JT2\. The initial particle velocity is 
assigned according to Eq. ( 133|) by setting t = and U = 1. Same as in [12], the 
initial configuration with 3600 particle is taken from previously stored particle 
positions (relaxed configuration). 

Figure [T] shows the calculated decay of the maximum velocity and the vor- 
ticity field at t = 1. It can be observed that the present solution is in quite 
good agreement with the exact and reference solutions fL2\. This is not un- 
expected since the only notable difference between the SPH-a model and the 
SPH discretization in [12] is the additional effective stress term in Eq. (!20|) . 
whose contribution is small since the flow is well resolved. Note that if the 
additional effective stress term is not applied, as expected, the numerical solu- 
tion (not shown here) has no noticeable difference with the reference solution. 
Note that, compared to the reference solution, the present solution is slightly 
less dissipative, suggested by the somewhat less errors in regions close to the 
centers of vortex cells. This behavior suggests that, as will also be shown in the 
next case, the overall effect of the additional effective stress term is decreasing 
dissipation. 
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4-2 Two dimensional decay turbulence 

We consider a two dimensional decay turbulence with high Reynolds number. 
The simulations are carried out with two resolutions, i.e. with 50 x 50 and 
100 X 100 particles. The particles were placed initially on a grid of squares in 
motion with velocities specified by a 8 x 8 array of Taylor-Green vortices 

u = cos(87rx) sin(87r?/), v = sin(87rx) cos(87r?/). (34) 

This initial condition mimics the experimental setup in which arrays of vortices 
are generated by applying electric field [23], or movable solid bars |15] . 

The energy spectrum at t = 2, 5 and 15, and the evolution of the normalized 
total enstrophy are given in Fig. |2l As shown in Fig. [2^, full turbulence has 
been developed at t = 2, and is characterized by a typical direct entropy cas- 
cade with —3 scaling of the energy spectrum [26] independent of the resolution 
of simulation. The spectrum at high wave numbers is that of white noise. This 
reflects the SPH method's limitation on resolving flow structured at small 
scales close to the smoothing length. Similar phenomenon is also observed in 
experiments, such as in [22j, where the measured spectrum is dominated by 
white noise at the high wave numbers beyond resolvability. As shown in Fig. 
Ub, the overall decay of enstrophy scales at about t^°-^ — independent of 
the resolution of simulation, which is in good agreement with previous direct 
numerical simulation [2] and theoretical prediction [28]. Note that if the addi- 
tional effective stress term in Eq. f l2U]) is not included in the model, as shown 
in Fig. [2]d, the enstrophy decay is considerably faster in the low resolution 
simulation. One important property of two-dimensional decay turbulence is 
the merging and pairing of vortices. As the velocity flelds at t = 2 and 15 
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shown in Fig. |3l the initially smaller vortices finally develop to the config- 
uration with a pair of big vortices comparable to square-size with opposite 
sign, which is known as the equilibrium state of two-dimensional decay tur- 
bulence [26]. This change is also refiected in the energy spectrum, as shown 
in Fig. m where the energy spectrum is considerably flatten. The probability 
density function (PDF) of velocity increments along the particle trajectory or 
acceleration at t = 2 and 15 is shown Fig. |H While small acceleration fluc- 
tuations follow a Gaussian distribution, large acceleration fluctuations show 
a distinct non-Gaussian. This PDF suggests the intermittency of the velocity 
fleld, which is well established for turbulence. Note that, while the turbulence 
is approaching the equilibrium state, other than disappear, the intensity of 
intermittency increases slightly. 

4-3 Two dimensional forced turbulence 

We consider a two dimensional forced turbulence with high Reynolds number. 
The simulation is carried out with 100 x 100 particles, which were placed 
initially on a grid of squares with zero velocity. The driving force is pairwise 
and produced by an inverse-viscous term 

= ^'™^e,, ■ VW,„ (35) 

where rj'^'^^ = 2.5 x 10^"^ and the velocity proflle v*'"" is obtained from the 
equation 

= cos(167rx+0) sin(167r?/+0), w"^™ = sin(167rx + 0) cos{16ny + (p) , (36) 

where the phase shift is produced with increments of a Wiener process. 
Again, this driving force mimics the experimental setup which drives an arrays 
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of vortices by applying electric field [23] . The difference is that zero mean flow 
is achieved strictly in the simulation. 

The energy spectrum at t = 2, 5 and 15, and the evolution of the total energy 
and total enstrophy (both are normalized with the values at t = 0.5) are 
given in Fig. |3 As shown in Fig. full turbulence has been developed at 
t = 2, in which depending on smaller or larger than the injection wave number, 
the energy spectrum show a direct entropy cascade with —3 scaling and an 
inverse energy cascade, which suggests transporting kinetic energy to larger 
length scales, with —5/3 scaling [26]. As there is energy injection, the kinetic 
energy and enstrophy increase at t^'^ before about t = 5. After this time, 
while the enstrophy increase is slightly slowed down to about the kinetic 
energy increase is slightly sped up to t^-^. This change is corresponding to 
the time, also see in Fig. [5^, from which the energy with the largest possible 
scale (square size) piles up. As shown in Fig. [5^, at t = 15 the turbulence 
reaches the so-called condensed state, where the energy at the largest scale 
is considerable large than that represented by a —5/3 inverse-cascade scaling 
and still increases with time. As shown in Fig. El while the velocity field at 
the early stage of forced turbulence resembles that of the decay turbulence, 
it is very different at the late time, when the forced turbulence reaches the 
condensed state and the decay turbulence the equilibrium state, though both 
are dominated with large-scale flow structures. In the condensed sate, different 
from the decay turbulence, there is no large vortex but large scale shearing 
flow. In addition, there are small vorticies corresponding the injection wave 
number, but no small scale vorticies found in the equilibrium sate of a decay 
turbulence, as shown in Fig. [3]d. The PDF of acceleration at t = 2, 5 and 15 is 
shown Fig. [7^. It can be observed that before the condensed state, the PDF is 
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in good agreement with that of decay turbulence before the equihbrium state. 
Compared to the decay turbulence in the equilibrium state, the intensity of the 
intermittency increases considerably in the condensed state, which is indicated 
by the narrowed probability for small fluctuation and extended probability for 
large fluctuation [26]. This is also reflected in the vorticity contour of Fig. [7b, 
where a number of spatially highly intermittent small pockets are presented 
with large positive and negative vorticity. 

5 Concluding remarks 

We have proposed a coarse-grained particle (SPH-cr) model for incompressible 
NS equation based on spatial filtering by utilizing SPH approximations. The 
SPH-(T particle model is similar to the LANS-a model and the SPH-e model 
but with the different additional effective stress term and approach to obtain 
the particle transport velocity. Since, numerically, this model resemble to the 
discretization of our previous developed incompressible SPH method, the same 
numerical method is applied. Numerical tests on two-dimensional moderate 
flow, and decay and forced turbulences are carried out. The results suggest that 
this model can be used for simulating incompressible turbulent flow problems 
to which SPH has been applied. 
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(a) t'lne (b) X 

Fig. 1. Taylor-Green problem simulated with 3600 particles: (a) decay of the max- 
imum velocity (reference: solution in [12j), (b) vorticity field (dashed line: solution 
in [12], dash-dot line: SPH-a) and exact solution (solid line). 




Fig. 2. Two dimensional decay turbulence: (a) energy spectrum at t = 2, 5 and 15. 
(b) evolution of total kinetic energy and enstrophy. The computations are performed 
on the periodic domain of [0, 0] x [1, 1] with 50 x 50 and 100 x 100 particles. 
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Fig. 5. Two dimensional forced turbulence: (a) energy spectrum at i = 2, 5 and 
15. (b) evolution of total kinetic energy and total enstrophy. The computations are 
performed on the periodic domain of [0, 0] x [1, 1] with 100 x 100 particles. 




X 



Fig. 6. Two dimensional forced turbulence: velocity field at {a) t = 2 and (b) t = 15. 
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Fig. 7. Two dimensional forced turbulence: (a) PDF of acceleration at i = 2, 5 and 
15; (b) 15 vorticity contours from —50 to 50 at t = 15. 
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